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Black hole-neutron star (BHNS) binaries are expected to be among the leading sources of grav- 
itational waves observable by ground-based detectors, and may be the progenitors of short-hard 
gamma ray bursts (SGRBs) as well. We discuss our new fully general relativistic calculations of 
merging BHNS binaries, which use high-accuracy, low-eccentricity, conformal thin-sandwich con- 
figurations as initial data. Our evolutions are performed using the moving puncture method and 
include a fully relativistic, high-resolution shock-capturing hydrodynamics treatment. Focusing on 
systems in which the neutron star is irrotational and the black hole is nonspinning with a 3:1 mass 
ratio, we investigate the inspiral, merger, and disk formation in the system. We find that the vast 
majority of material is promptly accreted and no more than 3% of the neutron star's rest mass is 
ejected into a tenuous, gravitationally bound disk. We find similar results for mass ratios of 2:1 
and 1:1, even when we reduce the NS compaction in the 2:1 mass ratio case. These ambient disks 
reach temperatures suitable for triggering SGRBs, but their masses may be too small to produce the 
required total energy output. We measure gravitational waveforms and compute the effective strain 
in frequency space, finding measurable differences between our waveforms and those produced by 
binary black hole mergers within the advanced LIGO band. These differences appear at frequencies 
corresponding to the emission that occurs when the NS is tidally disrupted and accreted by the 
black hole. The resulting information about the radius of the neutron star may be used to constrain 
the neutron star equation of state. 

PACS numbers: 04.25.D-,04.25.dk,04.30.-w 



I. INTRODUCTION 

Mergers of compact binaries, consisting cither of neu- 
tron stars (NS) or black holes (BH), are expected to 
be among the most promising sources of gravitational 
waves detectable by ground-based laser interferometers 
like LIGO Q, Q, VIRGO 0, i], GEO @, and TAMA 
[1, 0], as well as by the proposed space-based interfer- 
ometers LISA @ and DECIGO %. Theoretical mod- 
els indicate that a neutron star-neutron star (NSNS) 
[ToL [Til . [l2l IT3I [Til JlEI. ITHl or black hole-neutron star 
(BHNS) [ll Hilll Ha THj, S3 merger may result in 
a hot, massive disk around a BH, whose temperatures 
and densities could be sufficient to trigger a short-hard 
gamma-ray burst (SGRB). Indeed, SGRBs have been re- 
peatedly associated with galaxies with extremely low star 
formation rates (see [22| and references therein for a re- 
view), indicating that the source is likely to involve an 
evolved population, rather than main sequence stars. 

Modeling the inspiral, coalescence and merger of com- 
pact binaries requires fully general relativistic dynamical 
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simulations, and has been a long-standing goal of numer- 
ical relativity (see [23| for a review) . Historically, the first 
successful dynamical simulations of compact binaries in- 
volved NSNS binaries [jj, M, M, HI M, M, M, SI)- 
A breakthrough in the simulations of BHBH binaries oc- 
curred more recently [29|, [3(| HT| . Simulations of BHNS 
binaries, on the other hand, have so far lagged behind 
- perhaps because they combine the difficulties associ- 
ated with black hole singularities with the subtleties of 
relativistic hydrodynamics. To date the only fully self- 
consistent dynamical simulations of BHNS inspiral and 
coalescence are those of Shibata and Uryu pjj [2(| (here- 
after SU) and Shibata and Taniguchi 121] (hereafter ST). 

Over the past years, we have systematically developed 
the tools necessary to simulate the inspiral and merger 
of BHNS binaries, including the tidal disruption of the 
neutron stars and the potential formation of an accretion 
disk. As reviewed below, we have constructed quasiequi- 
librium initial data descri bing relativistic BHNS binaries 
in quasicircular orbits [32 . l33l . 1 34 , l35l . 1 3 61 ] , and performed 
preliminary relativistic dynamical simulations by assum- 
ing several simplifying approximations [Tt], Gil . We also 
demonstrated and tested how the numerical techniques 
adopted in many recent BHBH puncture simulations can 
be combined with relativistic hydrodynamics to evolve 
quasiequilibrium initial data that are provided only in 
the black hole exterior [3?], [H| • In this paper we report 
on our first fully self-consistent, dynamical simulations of 
BHNS binaries. 

For the initial data, we have adopted a hierarchical 
approach to construct quasiequilibrium models of BHNS 
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in quasicircular orbit. We began by making a number 
of simplifying assumptions, and have relaxed these as- 
sumptions step by step (|1 111, cf. [H, El 
for other BHNS initial data). Our current models, which 
we adopt as initial data for the dynamical simulations 
described in this paper, are solutions to Einstein's con- 
straint equations in the conformal thin-sandwich (CTS) 
decomposition. We model the neutron star as an irrota- 
tional r = 2 polytrope, and impose the black hole equi- 
librium boundary conditions of Cook and Pfeiffer 41.] on 
the black hole horizon, to approximate an irrotational 
BH. In our most recent paper [3q . we adopted the meth- 
ods of Caudill et.al. [43] to construct irrotational black 
holes more accurately, and found closer agreement with 
post-Newtonian results. These improved initial data will 
be incorporated into our next set of dynamical calcula- 
tions. 

For the purpose of comparison we point out that the 
dynamical simulations of SU and ST adopt initial data 
that are different from ours. Both approaches lead to 
valid solutions to Einstein's constraint equations, but the 
solutions may be physically distinct. Specifically, our ini- 
tial data use the CTS decomposition, which allows us 
to impose an approximate helical Killing vector on the 
spacetime and thereby set to zero several time deriva- 
tives of the field variables in a corotating frame. For 
example, we impose all conditions dt^fij — 0, where jij 
is the conformally related spatial metric, and these im- 
mediately yield a relation between the components of the 
extrinsic curvature and the shift vector [see Eq. (|27|) ]. By 
contrast, SU and ST do not impose these conditions, but 
instead use the conformal transverse-traceless (CTT) de- 
composition to obtain the extrinsic curvature on the ini- 
tial slice. However, they do employ the assumption of a 
helical Killing vector to construct a lapse and shift, (cf., 
[43I ] who use a similar approach for BHBHs) and these 
gauge quantities are used to solve the quasiequilibrium 
fluid equations for the neutron star. They are also used 
to compute the matter source terms appearing in the 
constraint equations. For example, SU and ST take the 
divergence of Eq. ([27)) to generate three equations for the 
shift. In addition, they model the BH as a "puncture" 
(see [HI, |45|, [4(| ) , whereas we excise the BH interior and 
impose the equilibrium boundary conditions on the ap- 
parent horizon to force the BH to be stationary, at least 
momentarily. Therefore, one might speculate that our 
CTS initial data may represent quasiequilibrium BHNSs 
in quasicircular orbit more faithfully than SU and ST's 
initial data. In this paper, we present evidence that the 
details of the initial data have a noticeable impact on the 
outcome of the merger - including the disk mass - which 
may explain some of the differences between the findings 
of SU and ST and ours. 

In most dynamical simulations of BHNS binaries to 
date, the self-gravity of the NS and/or the tidal gravity 
of the BH are treated in a Newtonian or post-Newtonian 
framework (see, e.g., 47, 48, 49, 50, 511; see also 1521 
who performed fully relativistic simulations of head-on 



collisions). In many calculations, especially those with 
an initial mass ratio q — Mbk/Mns ^5 3, significant disks 
are formed after the NS is disrupted, and for very stiff 
nuclear equations of state (EOSs), the core of the NS 
may survive the initial mass transfer episode and remain 
bound. These findings contrast with some semi-analytic 
relativistic arguments that suggest that it is very difficult 
to form disks with appreciable masses in the merger of 
BHNS binaries [H. 

Using our earlier initial data [12, , which assumed 
extreme mass ratios with g > 1, we performed simu- 
lations of BHNS merger in an approximate relativistic 
framework [l7], EH- In particular, we assumed that the 
spatial metric remains conformally flat throughout the 
evolution (see 0, HH). Though this approach only al- 
lows for crude estimates, we found that mergers of irro- 
tational BHNS binaries may lead to disks of masses up 
to 0.3 Mq, with sufficient heating to emit the neutrino 
fluxes that are required to launch a gamma-ray burst. In 
their fully relativistic BHNS simulations, SU later found 
disk masses in the range of 0.1 - 0.3 Mq for corotating 
NSs, and ST found smaller disk masses of 0.04 - 0.16 M 
for more realistic irrotational NSs. 

In preparation for our fully relativistic dynamical sim- 
ulations of BHNS merger, we demonstrated in [3?], [38j] 
that the moving puncture method (see [3(1 HH as well 
as numerous later publications), which has proven ex- 
tremely useful for BHBH simulations, can be adopted for 
BHNS simulations and conformal thin-sandwich initial 
data. Two conceptional issues were addressed, namely 
the inclusion of relativistic hydrodynamics into these sim- 
ulations, and the fact that moving puncture simulations 
require initial data everywhere, while our conformal thin- 
sandwich initial data excise the black hole interior and 
hence provide data only in the black hole exterior. 

In contrast to the original dynamical puncture simula- 
tions [H, H3] , in which the puncture was forced to remain 
at a fixed coordinate location, the "moving puncture" ap- 
proach allows the punctures to move freely through the 
computational grid. This method is typically used in 
the context of the BSSN formulation [H, [59|, coupled 
to a "Gamma-driving" shift [6(| and a "1+log" slicing 
condition 6l|. Geometrical arguments show that, with 
this slicing condition, dynamical simulations approach 
limit surfaces of finite areal radius around black hole 
sing ularities, but never reach the singularity itself (see 
[621163. HH, [65J ) . These findings provide insight into why 
moving puncture simulations can possibly be successful, 
and also suggest that it may be possible to incorporate 
relativistic hydrodynamics into these simulations. Since 
the simulations only cover regular regions of the space- 
time, the hydrodynamic flow never encounters any black 
hole singularities. In [3?| we demonstrated that with 
only very minor modifications, our high-resolution shock- 
capturing (HRSC) relativistic hydrodynamics algorithm 
(see (66|) can indeed be used together with the moving 
puncture method to model accretion onto black holes. 

Part of the appeal of the moving puncture approach 



3 



stems from the fact that it does not require an excision of 
the black hole interior. Accordingly, this method requires 
initial data everywhere, both in the BH's exterior and in- 
terior. Most dynamical moving puncture simulations to 
date have therefore adopted the puncture method also 
in the construction of the initial data (see e.g., [46L 
As discussed above, CTS initial data are generally be- 
lieved to be very good approximations of true quasiequi- 
librium states, but solving the CTS equations usually 
involves excising the BH interior, so that the resulting 
data exist only in the BH exterior. By definition, no 
physical information can propagate from the BH interior 
to the exterior, and we have recently demonstrated that 
even unphysical, constraint- violating noise ("junk") does 
not leave the BH ((37l. l38l |. see also (68[) in numerical 
evolutions employing the BSSN formulation and moving 
puncture gauge conditions (modulo certain smoothness 
restrictions on the junk data near the horizon). Thus, 
the BH interiors in CTS initial data can be filled with 
"junk" without affecting the external spacctime, enabling 
us to evolve our quasiequilibrium BHNS initial data via 
the moving puncture formalism. 

With all the aforementioned pieces in place, we now 
report our first fully self-consistent, relativistic dynam- 
ical simulations of BHNS binaries. We are particularly 
interested in binaries that may potentially lead to a siz- 
able accretion disk, so we focus on binaries in which the 
neutron star is tidally disrupted just before reaching the 
innermost stable circular orbit (ISCO) and plunging into 
the black hole. The binary separation c? t id at which the 
neutron star will be tidally disrupted may be estimated 
from the following crude Newtonian argument. Equating 
the tidal force exerted by the BH on a test mass at the 
NS's surface with the gravitational force exerted by the 
neutron star on this test mass, we find 



tid 



Mi 



-2/3 C - 



BH 



(1) 



where C = Mns/Rns is the neutron star compaction. 
Given typical neutron star compactions of C ~ 0.2, small 
but reasonable values of q arc required for dtid to be larger 
than the ISCO separation of about disco ~ 6A/bh- Our 
more careful analysis ([III; see Fig. 15) shows that for 
r = 2 polytropes we need q < 4.25 . Given typical 
neutron star masses (Mjvjs ~ 1.5M Q ), this means that 
we can expect the formation of an accretion disk only for 
low-mass black holes. 

How often such binaries merge in the observable uni- 
verse is still an open question. The uncertainties arise 
from some aspects of population synthesis calculations 
that are only partially understood. In particular, en- 
velope ejection efficiency during the common envelope 
phase seems to be a crucial factor in forming low-mass 
BHs during binary stellar evolution. For example, if one 
assumes efficient ejection and a large maximum NS mass, 
the primary NS will generally accrete insufficient mass to 
induce collapse to a BH, and one ends up with a large 
number of NSNS binaries. For inefficient ejection and 



a smaller maximum NS mass, it is relatively easy for 
the NS to accrete sufficient material to form a BH with 
a mass only slightly larger than a NS. Although some 
previous population synthesis calculations working un- 
der the latter assumption found a nearly flat spectrum 
of binary mass ratios spanning the range q = 1.5 — 10 
[6^ |. a more recent calculation that assumes highly effi- 
cient envelope ejection yields typical binary mass ratios 
q = 6 — 10 70]. The latter scenario would predict that 
NSs undergoing tidal breakup prior to reaching the ISCO 
are rare, as are any resulting SGRBs from these systems. 
The overall rate estimates for BHNS mergers observable 
by an advanced LIGO detector typically fall in the range 
TZ ~ l-lOOyr" 1 [n|- 



These issues noted, we begin our investigation of BHNS 
binary merger and coalescence in full general relativity. 
This paper is the first in a sequence of papers which will 
thoroughly explore the effect of various binary parame- 
ters on the tidal disruption, disk formation, the poten- 
tial for launching a GRB, and the corresponding grav- 
itational wave signals. In this paper we will focus on 
irrotational binaries with mass ratios q < 3. As an addi- 
tional word of caution, we point out that our results are 
fundamentally limited by uncertainties about the true 
nuclear EOS, both in the cold initial state as well as the 
later shock-heated hot phase. Disk masses may depend 
sensitively on the structure of the NS, especially the low- 
density outer regions, so all BHNS merger results should 
be viewed in light of this caveat. In particular, the like- 
lihood of BHNS mergers as SGRB progenitors may be 
difficult to determine conclusively until this issue is re- 
solved. 

This paper is organized as follows. In Sees. |TT] and 
IIII1 we summarize the basic equations and their specific 
implementation in our general relativistic hydrodynam- 
ics scheme, along with a discussion of initial data, gauge 
conditions, matter evolution, and diagnostics. In Sec. HVI 
we discuss the results of our BHNS merger simulations, 
and how they depend on both physical as well as compu- 
tational parameter choices. We conclude in Sec. [V] with 
a discussion of our findings, and comment on future di- 
rections. 



II. BASIC EQUATIONS 

In this section we list the full set of evolution equa- 
tions integrated by our numerical code. Field, coordi- 
nate, and hydrodynamic evolution equations are sum- 
marized in Sees. Ill Al III B[ and III CI respectively. 



A. Field Evolution: The BSSN Equations 

Assuming geometrized units in which G = c = 1, we 
write the spacetime metric in the standard 3+1 form 



ds 2 



a 2 dt 2 + jij(dx l + p'd^idx 3 + p>dt), (2) 
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where a, and 7^ are the lapse, shift, and spatial 3- 
metric, respectively. The extrinsic curvature K t j is de- 
fined by 



{d t - £phij = -2aK, 



(3) 



Here Cp is the Lie derivative with respect to j3 l . 

In the BSSN formalism, we define the conformally re- 
lated metric 7^, the conformal exponent 0, the trace of 
the extrinsic curvature K, the conformal traceless extrin- 
sic curvature A i2 - , and the conformal connection functions 
as follows 



4> 


= ^ln[det( 7lJ )], 


(4) 




= e-'Hj, 


(5) 


K 


= 7ijK ij , 


(6) 






(7) 


r 


= ~f j J, 


(8) 



where j denotes the partial derivative: j lJ j = d^ % K We 
use the same field evolution equations as Eqs. (11)— (15) 
of O: 



(d t - Cp)lij = -2aAij, 
(dt-CpW = -~aK, 

[d t -C p )K = -^D^a+^aK 2 

+aA lj A i: > +47ra(p+ S), 



(9) 
(10) 

(11) 



{d t - C p )Aij = e-^i-DtDja + a{R %3 - 8^)) TF 
+a(KA t3 - 2A il A l j ), (12) 

and 

d t r = dj(2aA ij + £pf j ) 



(13) 



-2a [ -f 3 Kj - 6A ll 4>j - r jk A 3k + H-irf^Sj 



This condition [37] 
a fisheye radius f, 
pression is closely 
of shift evolution 
that setting 77 ~ 
evolutions, where 
as defined in Eq. 
except run B, for 



I is similar to that in [3l| , but allows for 
discussed in Sec. IIIIAl below. This ex- 
related to the "Gamma-driver" family 
equations. We have found empirically 
0.5 /M yields well-behaved coordinate 
M is the ADM mass of the system, 
(f3"T)) . This value is chosen for all runs 
which we use 77 = 0.413/M. 



C. Hydrodynamic Equations 

The matter source terms are defined as 

p = n a nf3T a(i , 
Si = ~l la n & T^ , 
Sij — "fi a "fjpT al3 , 



(17) 



where T Q/3 = (p + p e + P)u a u p + Pg a0 is the stress- 
energy tensor for a perfect fluid, po, e, P and u a are the 
fluid's rest-mass density, specific internal energy, pres- 
sure, and 4-velocity, respectively, and n a = (— a, 0,0, 0) 
is the future-directed unit normal to the time slice. 
We evolve the "conserved hydrodynamic variables" , 
Si and f, defined as follows 

p* = -n^pou* 1 = ay/^pou , (18) 
Si = y/jSi = ay/jT°i = p*hui, (19) 
t = y /jn li n v T>"' -p. =a 2 ^T m -p*. (20) 

The evolution equations for these variables are given by 
Eqs. (34), (36), and (38) of [H, 



d t p* + d j (p tt v j ) = 0, 



(21) 



dtSi + djia^fT^) = X -a^T^d l9a ^(22) 
dtr + d^VlT 01 - P*v l ) = s, (23) 
where 7 = det(7y) = e 12c ^, and the energy source term s 



B. Gauge Equations 

As in most moving puncture calculations, we use an 
advective "1+log" slicing condition for the lapse 



d t a — f3 l diCt = 2aK, 



(14) 



and a second-order "non-shifting-shift" (in the language 
of |li) 



dtP = -B* 



dr 



dtB % = ( — ) d t T l - V B l 



(15) 
(16) 



s = -a^T^V v n^ 

= a/y [{T m pf3 J + 2T 0i (3 j + T i] )K l3 

-(T m p l + T 0l )d l a}. (24) 

Here v l = u l /u° is the fluid's 3-velocity. 

To complete the system of equations, we specify an 
EOS. Our code is capable of handling EOSs of the form 
P = P(p ,e). In this paper, we employ the standard 
T-lawEOS 



P = (T - l)poe 
with r = 2 to model the NS matter. 



(25) 
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III. NUMERICAL METHODS 

The code we use is very similar to that described in 
[37l [38j . We do not consider magnetic fields in this 
paper, so the magnetic field sector is disabled in these 
calculations. The equations of general relativistic (GR) 
hydrodynamics are handled by a HRSC technique [661 ] 
that employs the monotonized central (MC) reconstruc- 
tion scheme coupled to the HLL (Harten, Lax, and 
van Leer) approximate Riemann solver 17511 . The metric 
is evolved via the BSSN formalism (5f| [5!| as described 
in [72| . but with fourth-order accurate spatial differenc- 
ing and upwinding on the shift advection terms. Our 
code is based on the Cactus parallelization framework 
[76j , in which our second-order Iterated Crank- Nicholson 
time-stepping is managed by the MoL, or method of lines, 
thorn. 

For completeness, we provide below a brief overview 
of our grid setup and initial data (Sec. IIII A[) ; a dis- 
cussion of field, gauge (Sec. IIII Bp . and hydrodynamic 
( Sec. IIII CI evolution techniques; a description of how we 
apply boundary conditions (Sec. IIII Dj) : a summary of 
diagnostic techniques we use to both validate our numer- 
ical results and examine our spacetimes fSec. IIII E|) ; and 
finally a description of the technique we use to measure 
gravitational wave (GW) emission fSec. IMF]) . 



A. Grid Setup and Initial data 

We use a "fisheye" coordinate system [77j to expand 
the physical extent of our numerical grid while main- 
taining high resolution in the strong-field region. Signif- 
icantly lower resolution is maintained in the wavezone, 
but it is set so that a gravitational wavelength is resolved 
by at least 12 gridpoints. To set up a fisheye coordinate 
grid, we define the "physical" radius r in terms of a fish- 
eye radius f according to 

r - a f + ( a »-i - a i) s i , ( cosh((r + _R i )/s i ) \ 
r-a n r ^ 2 tanh( J R i /si) n v,cosh((f - R l )/s l )) ' 

(26) 

Here di sets the magnitude of the i'th fisheye transition, 
Si determines the width of the transition, and Ri specifies 
the center of the transition. In physical coordinates, the 
grid spacing smoothly transforms from Ax « a^iAx to 
Aa; ss a.iAx, over a set of fisheye coordinates spanning 
radii Ri — Si < f < R4 + s^. For convenience, we always 
set ao = 1, so that our coordinate grid spacing Aa; rep- 
resents the physical grid spacing in the central region of 
our numerical grid. In this paper we use only one tran- 
sition zone (n = 1), with a\ — 8. We have listed other 
relevant grid parameters in Table ITT1 

We defer to Appendix A of [37| for transformation laws 
pertaining to all field and hydrodynamic quantities un- 
der a fisheye transformation. All of our calculations are 
performed assuming equatorial symmetry, on numerical 



grids of the form 2N x 2N x N, with N ranging from 166 
to 305. 

All initial data we evolve in this paper were generated 
by [351 ] - To map these spectral configurations onto our 
non-spectral simulation grid, we first construct our nu- 
merical grid and record the positions of each point in 
physical coordinates. Then we evaluate the field and 
hydrodynamic quantities based on their spectral coeffi- 
cients. Next, we transform the vector and tensor quan- 
tities into fisheye coordinates via transformations found 
in Appendix A of 37]. Finally, the excised BH region 
is filled with constraint-violating initial data, using the 
"smooth junk" technique we developed and validated in 

M 

The assumptions under which our initial data are 
constructed differ from those of SU and ST. We solve 
Einstein's constraint equations in the conformal thin- 
sandwich (CTS) formalism, which allows us to impose 
an approximate helical Killing vector by setting the time 
derivatives of the conformally related metric to zero. As 
a result, our extrinsic curvature is always related to the 
shift vector that appears in the solution through Eq. (4) 

of m, 

i« = ^ + - V fc /3 fc ^ . (27) 

These quasiequilibrium initial data excise the black hole 
interior, allowing us to impose equilibrium boundary con- 
ditions on the excision surface. By contrast, SU and ST 
adopt the CTT decomposition to obtain the extrinsic 
curvature on the initial slice, but employ the assump- 
tion of a helical Killing vector to construct a lapse and 
shift. Also, they model the black hole as a puncture (see 
[43l HI, HI, |46|]). Both sets of initial data are solutions 
to Einstein's constraint equations, but they may differ in 
both the amount of spurious gravitational wave content 
and the degree of orbital eccentricity. 

B. Metric Evolution and Gauge 

We apply two methods that have improved the stabil- 
ity and accuracy of our field and gauge evolution when 
evolving BH spacetimes. 

First, we use fourth-order finite differencing schemes to 
calculate spatial derivatives in the field/gauge evolution 
sectors. Also, for any terms of the form j3 l di . . ., which 
arise in both the Lie derivative terms and in the lapse evo- 
lution, we use fourth order upwind differencing stencils 
instead of the standard centered fourth-order stencils (see 
Eqs. (2.5)-(2.6) and (2.2)-(2.4) of [z|, respectively). We 
note for completeness that our mixed second-derivative 
stencil is slightly different that given by Eq.(2.4) of [78| . 
but remains fourth-order convergent, 

d X yFi,j,k = 4g^ x i F i-2,j+2,k + F l+2 ,j-2,k (28) 
— Fi+2,j+2,k — Fi_2,j-2,k 
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+16(i 7 i+i,j+i,fe + -Fi-i,j-i,fc 
—Fi-lj+l,k ~ Fi+l,j-l,k)] ■ 

In addition, wc enforce the conditions 7 = det(7y) 
I and A = tr (Ay) = at every timestep, using the 
substitutions 



Hi 



An 



(29) 
(30) 



our grid, with a density floor set to 10~ 10 of the maxi- 
mum density on our grid at t = 0. The initial atmo- 
spheric pressure P atm is set to the cold polytropic value 
-Patm = K Patm; where k is the polytropic constant at 
t = 0. Throughout the evolution, wc impose limits on 
the atmospheric pressure to prevent spurious heating and 
negative values of the internal energy e. 



D. Boundary Conditions 



as is commonly done in numerical relativity codes. We 
do not, however, enforce the Hamiltonian, momentum, 
or Gamma constraints, 



= H = f J D i D j e' p - 



-R 

,50 



(31) 



= M* 

o = g l 



= Dj(e 64, A ji ) - -e 6ct >D l K - Sne^S*, (32) 



r + 7 y 



(33) 



so these serve as an independent check on the validity 
of our code. We do not add a Hamiltonian constraint 
damping term to the right hand side (RHS) of the <j> 
evolution equation, but we add damping terms to the 
RHS of the T'\ jij, and A^ BSSN evolution equations, 
following the prescription defined in [79[ . These terms are 
zero analytically, and serve only to stabilize evolutions. 



We apply Sommerfeld outgoing wave boundary condi- 
tions to the entire set of field and gauge variables f 



f(r,t) 



-At 



f (r -Ar,t- AT) 



(36) 



on the outer boundary of our numerical grid. Here AT is 
the timestep and Ar = ae~ 2 ^ AT , where radii are evalu- 
ated in physical (as opposed to fisheye) coordinates. To 
enforce these boundary conditions, we first transform /3 Z , 
4>i Aij, and gij into physical coordinates, apply the Som- 
merfeld condition as specified above, and then transform 
back to fisheye coordinates. We do not transform T l , as- 
suming that its value propagates outward independent of 
the radius. 

Since we adopt the moving puncture method (as op- 
posed to black hole excision) there are no interior bound- 
aries or boundary conditions. 



C. Hydrodynamic evolution 

The hydrodynamics equations are calculated using the 
HRSC scheme described by 66]. To recover the "prim- 
itive variables" po, P, and v i from the conserved set 
p*, f, and Si, we perform the inversion as specified by 
Eqs. (57)— (62) of [66| . As in (3?j], our inversion algorithm 
occasionally finds unphysical sets of conserved variables 
at points immediately adjacent to the puncture and in 
our atmosphere, which do not allow for solutions of the 
primitive variables. As in that paper, we enforce the 
following two conditions, which are both necessary and 
sufficient to allow for a well-defined inversion everywhere, 
and result in smooth hydrodynamic variable profiles in 
the BH interior after the puncture has passed through a 
set of grid points: 

|5| 2 = 7 u '5 l S' J < f(f + 2p„), (34) 
r > 0. (35) 

When these conditions are not met we rescale S l so that 
its new magnigude is IS^ 2 = 0.98f(f + 2p*), and set f = 
fo-maxj where TQ- max is the maximum value of f 
present in our initial data. 

To stabilize our hydrodynamic scheme in regions where 
there is no matter, we maintain a tenuous atmosphere on 



E. Diagnostics 

To validate our calculations, we compute surface inte- 
grals for the system ADM mass M, linear momentum Pj, 
and angular momentum J.;, given by 

M = hf (g^-f^^)^, (37) 

p i = ^f(Ki-SiK)d^, (38) 

J i = h Cijk f xj ( K k-^k K ) d ^m, (39) 

where ip = and cE^ = (x z /r)ip e r 2 sin 9dddip for a 
spherical surface at fixed radius. Note that the above 
expressions are valid only if the spatial 3-metric 7^ ap- 
proaches the Minkowski metric rjij at large r. Hence we 
need to transform the variables from fisheye to physical 
coordinates before integrating. In addition, we monitor 
the following normalized expressions for the Hamiltonian 
and momentum constraints 

||W|| = / (\H\/N w )dV (40) 

\\M l \\ = f (\M l \/N MC )dV, (41) 
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where 



D l Di%jj 



^A 13 aA + (^k 2 

8 3 / V 12 




(42) 



N MC = [J2 



inS 1 ) 



< >2 



(43) 



+ 



U- 6 Dj{ij; 6 A ij ) 



1/2 



The apparent horizon (AH) of the BH is computed 
using the ahf inderdirect Cactus thorn (80j |. This thorn 
outputs the BH irreducible mass, which is related to the 
AH area A as follows: 



Mir 



y/A/lGir. 



(44) 



With the AH surface computed by ahf inderdirect, we 
can evaluate diagnostic integrals in a region interior or 
exterior to the AH. 



F. Gravitational radiation 

To measure the gravitational wave (GW) emission 
from our binaries, we use both gauge-invariant theory 
based on perturbations of a background Schwarzschild 
spacetime derived by Zerilli [§l[ and Moncrief [82[ , here- 
after referred to as the "Z-M" formalism, as well as 
a technique that makes use of the Newman-Penrose 
Weyl scalar ip4. In the Z-M formulation, deformations 
of the spatial metric are viewed as perturbations on a 
Schwarzschild background at large radii. We decom- 
pose these perturbations into gauge-independent even 
and odd-parity modes, denoted ^ l J£ en and Q 1 ^ in the 
notation of [831 ] . whose derivation is outlined below. For 
convenience, we define the time- integral of the odd-parity 
mode amplitude 



odd 



(45) 



In terms of these expressions, the complex gravitational 
wave strain H = h + — ih x is given by Eq. (4.34) of (83|, 



H 




+ i& 



Ira 
odd 



! iii 



(46) 



where _2^' m is the s = —2 spin-weighted spherical har- 
monic. In the notation of Shibata and collaborators (see, 
e.g., [HI and earlier papers), Rf m and Rf rn are related to 
our adopted notation as follows: 



V2r\ (1-2)! 



odd' 



(48) 



In addition, we use the PsiKadelia thorn to compute 
the complex Weyl scalar ^4, which depends on the spa- 
tial metric and extrinsic curvature. The wave strains are 
given in terms of ip4 by Eq. (3.4) of 83], 



H 



tpidt'dt". 



(49) 



We also decompose ip4 into s = — 2 spin- weighted spher- 
ical harmonic modes. 

In terms of H, the radiated energy, linear momentum, 
and z-component of the angular momentum are calcu- 
lated from Eqs. (2.8), (2.11), and (2.13) of H 



dE 
~dt 
dP i 
~~dT 

dJ z 
dt 



lim — I \H\ 2 dfl, 

r-+oo 167T / 



lim 

i — >oo 167T 



H\ 2 dtt, 



lim 



32tt 



H*d <l> Hd£l. 



(50) 
(51) 
(52) 



In practice, radiative losses are computed as a sum over 
modes (up to and including 1 = 4), following expressions 
equivalent to Eqs. (4.41), (4.43), and (4.47) of H for 
the Z-M case, and Eqs. (3.6), (3.8), (3.14), and (3.24) for 
the ip4 case. 

To compute the gravitational wave energy spectra, we 
use the same techniques as [14|, determining the energy 
loss per unit frequency from their Eq. (8), 



dE 



(Z-2)! 1 



/;;/ | 



where \^i. 



(53) 



and ^ 



| e 2vi f t '$>(t)dt. We also define the effective gravitational 
wave amplitude, taking the z — > limit of Eq. (5.1) from 



hes(f) 




(54) 



Notice that this expression differs by a proportionality 
constant from that found in [8(| . This is due to different 
assumptions regarding the alignment of the binary and 
viewing angle relative to the detectors. The constant 
here is specified to reflect the RMS average of signal am- 
plitudes over all possible orientations of merging binaries 
and the detector. 



IV. NUMERICAL CALCULATIONS OF BHNS 
MERGER 



nE _ 



1 



'(1 + 2) 



V2r V (I - 2) 



| even i 



(47) 



To begin our survey of the BHNS merger parameter 
space, we focus on those configurations likely to undergo 
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mass loss and disk formation, i.e., systems with binary 
mass ratios comparable to unity. To ensure that our re- 
sults are physically relevant, we consider cases with ap- 
propriately large NS compactions, noting that lowering 
the compactions would certainly increase the mass of a 
potential disk. 

All our calculations begin with initial configurations 
taken from [35[, which consist of irrotational NSs orbit- 
ing approximately nonspinning BHs. We investigate two 
different cases for the NS compaction, distinguished by 
nondimensional mass Mq = Mq/k 1 / 2 = 0.15 and 0.1, 
where M is the rest (baryon) mass of the NS, and k 
is the polytropic constant. We refer to these as the high 
and low-compaction NS cases. These configurations yield 
compactions C = M^s/Rns of 0.145 and 0.0879, respec- 
tively, where Mns is the ADM mass and i?NS is the radius 
of the NS in isolation. Setting k fixes the NS mass for 
a given compaction and polytropic index [§7|- Choos- 
ing the high (low) compaction NS to have a rest mass of 
1.4M , we find the ADM mass for the isolated NS to be 
1.30Af Q (1.34M©), with an isotropic radius of 11.24 km 
(20.46 km), and circumferential (Schwarzschild) radius of 
13.24 km and (22.49 km). 

We consider data from four different independent se- 
quences in [35] , consisting of the high-compaction NS in 
binaries with mass ratios q = Mj rr /Mjvs = 3, 2, and 
1, as described in their Table IV, along with the low- 
compaction NS for a binary mass ratio q = 2, from their 
Table V. Here, Mi rr is the BH irreducible mass. These se- 
quences are denoted "A"-"D", respectively, and are sum- 
marized in Table [U 

For sequence A, we consider three separate initial bi- 
nary separations (denoted Al, A2, and A3) to determine 
how the initial separation affects the results of the simu- 
lations. For case Al, we performed a series of simulations 
designed to study the effect of the numerical grid parame- 
ters on our results. In cases Al-lo, Al-med, and Al-hi, we 
varied only the grid resolution, keeping the outer bound- 
ary at a fixed location. Case Al-farbc has the same grid 
resolution as Al-med, but the outer boundary is farther 
away. 

Performing multiple runs with the same quasiequilib- 
rium binary parameters (sequence A) enables us to gauge 
the resources necessary to perform accurate evolutions 
and validate our physical results. For example, we may 
compare the location of the tidal disruption point against 
quasiequilibrium estimates of [36j. We also investigate 
how the disk mass varies with regard to our grid parame- 
ters, noting that this had a significant effect on the results 
found in [19l |20| . Finally, performing multiple runs al- 
lows us to check the robustness of the gravitational wave 
signals. 

For convenience, we identify the approximate time at 
which the low-density regions of the deformed NS first 
fall into the BH horizon as the moment of "first contact" , 
£fc- 



A. BHNS binaries with mass ratio 3:1 

To begin our discussion, we consider results from evo- 
lutions of sequence A, the high compaction NS case with 
M = 0.15 and q = 3. 

In Fig. [TJ we plot density contours with overlaid 3- 
velocity vectors in the equatorial plane for our large 
initial separation case (A3). The orbital direction is 
counter-clockwise in these snapshots. The upper left 
panel shows the initial configuration. We see the onset 
of accretion after approximately 1.75 orbits (t « 215M), 
with matter flowing in a narrow stream through the in- 
ner Lagrange point and into the BH. The accretion flow 
then accelerates as the NS is consumed by the BH. Later, 
at t « 290M, we see the beginnings of mass loss out- 
ward into a disk through the outer Lagrange point, but 
the mass stream is quite tenuous. At the end of our 
simulation (t — 418. 7M), no more than 3% of the total 
rest mass of the NS exists outside the AH. The remain- 
ing matter is gravitationally bound to the BH. However, 
not all of this matter will form what would typically be 
referred to as a disk, i.e., a quasi-stationary torus that 
evolves on secular rather than dynamical timescales. In- 
stead, some of the exterior mass at the end of our calcu- 
lations will be accreted on relatively short timescales, un- 
til the remaining fraction achieves rotationally-supported 
quasi-equilibrium. Thus, these estimates should be taken 
as upper limits on true "disk masses" for a particular set 
of initial data and grid parameters. 

In case A of [l|| , SU perform a fully GR simulation of 
a BHNS merger with a synchronized NS and a mass ratio 
of q = 2.47. They find the fraction of the initial rest mass 
of the NS in the final disk to be « 19 - 28% of the initial 
rest mass of the NS. We find this disk mass fraction to 
be < 3% for our q — 3, irrotational NS model A3, which 
suggests that more slowly spinning NSs feed less mass 
into a disk. In case A of 2l|, ST simulate an irrotational, 
q = 3.06 BHNS binary and obtain a disk with a rest 
mass that is » 6.6% of the initial NS rest mass. While 
this result is consistent with our observation that lower 
NS spin suppresses the mass of the disk, ST still find a 
significantly larger disk than we do. 

To check the validity of our simulations, we monitor 
the normalized Hamiltonian and momentum constraint 
violations using Eqs. (|40[) and (|4"Tj) . respectively. We 
show the results from sequence A in Fig. [51 In all cases, 
the Hamiltonian and momentum constraint violations are 
~ 3% throughout the evolution. We find that at early 
times, the Hamiltonian constraint violation decreases as 
resolution is increased, as expected. However, at late 
times the Hamiltonian constraint violation steadily in- 
creases, and its magnitude becomes insensitive to resolu- 
tion (see lines corresponding to cases Al-hi, Al-med and 
Al-low in Fig. suggesting that the error is dominated 
by reflection from the outer boundaries. Indeed, we find 
that the Hamiltonian constraint violation is significantly 
smaller when we move the outer boundaries outward (see 
case Al-farbc in Fig. [2]). 



9 



TABLE I: Summary of our initial configurations, which are taken from the results of [35l ]. Here, the mass ratio q = M\ tt /Mns, 
Mo = M /k 1/2 is the nondimensional rest (baryon) mass, M is the total ADM mass for the binary system (assuming that 
Mo = IAMq), clq the initial binary coordinate separation, J the initial angular momentum of the system, f2 the orbital 
frequ ency , and tFc the time of first contact. Finally, MSIcts is the frequency of tidal disruption, as derived from our initial 
data [351 ] . and MQ num is the frequency at which our numerically derived waveform spectrum deviates from the restricted 
post-Newtonian value by 25% (see Section TlV C|l . 



Case 


<7 


M 


M/Mq 


ao/M J/M 2 


Mfl 


tpc/M M^cts 




Al 


3.0 


0.15 


5.15 


5.41 0.629 


0.0628 


50 


0.0728 




A2 


3.0 


0.15 


5.16 


6.49 0.663 


0.0435 


105 


0.0728 




A3 


3.0 


0.15 


5.17 


8.81 0.698 


0.0329 


215 


0.0728 


0.0789 


B 


2.0 


0.15 


3.86 


7.17 0.774 


0.0499 


60 


0.0550 


0.0656 


C 


1.0 


0.15 


2.57 


8.61 0.936 


0.0343 


91 


0.0382 


0.0509 


D 


2.0 


0.1 


3.98 


11.56 0.909 


0.0228 


155 


0.0255 


0.0395 



TABLE II: Summary of the grid configurations and fisheye parameters used for our runs. 



Run 


Ax int /M Ax cxt /M 


Grid Size 


r out /M Ri/M 


si/M 


Al-hi 


0.05 


0.4 


540 2 


x 270 


39.5 


9.7 


1.6 


Al-med 


0.061 


0.49 


440 2 


x 220 


39.5 


9.7 


1.6 


Al-lo 


0.081 


0.65 


332 2 


x 166 


39.5 


9.7 


1.6 


Al-farbc 


0.061 


0.49 


604 2 


x 302 


79.6 


9.7 


1.6 


A2 


0.061 


0.49 


532 2 


x 266 


39.7 


12.9 


1.6 


A3 


0.061 


0.49 


516 2 


x 258 


40.0 


12.3 


1.6 


B 


0.081 


0.65 


404 2 


x 202 


50.2 


11.6 


2.2 


C 


0.05 


0.4 


610 2 


x 305 


40.0 


11.7 


1.6 


D 


0.061 


0.49 


566 2 


x 283 


40.1 


15.8 


1.6 



Next, we investigate how our results depend on both 
the initial binary separation and the numerical reso- 
lution. In Fig. [31 we plot the fraction of the rest 
(baryon) mass inside the apparent horizon, /j n = Mo(r < 
rAn)/Mo, as a function of time for sequence A. We find 
that fi n depends only very weakly on the resolution or 
the location of the outer boundaries. We do, however, 
find stronger dependence on the initial binary separation 
(column ao /M in Table HJ . This dependence may be due 
to the zero radial infall speed in our initial data, which re- 
sults in a slightly eccentric orbit that increases the radial 
infall speed as the binary approaches the ISCO, affect- 
ing the tidal disruption and the disk formation (compare 
[H [H, H [§3). This effect may be compounded by 
the fact that our initial configuration is very close to the 
ISCO. Alternatively, this dependence may be caused by 
growing numerical error during the simulation. Despite 
these small uncertainties, we find /; n ^ 97% at the end 
of our simulations in all our cases, suggesting that there 
is no appreciable disk left behind after merger. 

Figure [4] shows the evolution of the BH irreducible 
mass, Mi rr , for the sequence Al simulations. We see that 
Mi rr increases as the NS matter is accreted, as expected. 
At late times, M; rr approaches an asymptotic value when 
most of the matter has fallen into the BH. However, re- 
flection from the outer boundaries causes M- m to slowly 
increase at very late times for runs Al-hi, Al-med and 
Al-low. This spurious effect is significantly reduced when 



the outer boundaries are moved to a larger radius (run 
Al-farbc), and M irr — > 0.954M at late times. 



B. The effect of the binary mass ratio and NS 
compaction 

Fig-Eldemonstrates the variation in dynamics due to a 
change in binary mass ratio, keeping the NS companion 
fixed. All three cases evolve approximately 0.5 orbits be- 
fore the tidally disrupted NS first touches the apparent 
horizon (i.e, the time of "first contact" - the top 3 plots 
of Fig. [5]). At this point in time, the general morphology 
of the systems are the same, regardless of mass ratio: a 
funnel-shaped NS, with matter flowing through the nar- 
row end of the funnel into the BH. After the time of "first 
contact" {t 1 = t — tpc > 0) , the dynamics of the system 
depend most sensitively on the initial mass, and hence 
size, of the black hole in comparison to the size of the 
NS. Defining x f™» o1 as the (coordinate) angular extent 
of the accretion funnel (with the density cutoff as defined 
in Fig.Q} around the AH on the equatorial plane, we find 
that as X/ 1 " 1161 increases to 180°, the accretion rate slows. 

For example, in the A cases (q = 3), the BH is suffi- 
ciently large so that X^ nncl < 180° throughout much of 
the simulation. We see that the NS accretion rate does 
not slow down until ss 90% of the rest mass falls into the 
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FIG. 1: Snapshots of density and velocity profiles at selected times for run A3, with binary mass ratio q = 3. The contours 
represent the density in the orbital plane, plotted logarithmically with four contours per decade, with greyscaling added for 
clarity. Arrows represent the velocity field in the orbital plane. The minimum contour value in each frame is Kpo(min)= 10~ 4 , 
or p (min)= 7 x 10 11 (1.4M /M o ) 2 g cm -3 . The maximum initial NS density is npo = 0.13. We specify the black hole AH 
interior in each snapshot with a filled black circle. In cgs units, the total ADM mass for this case is M = 3 x 1O" 5 (M O /1.4M ) 
s= 8(M o /1.4M )km. 



BH (Fig. ©. After the densest part of the NS falls into 
the BH (t' > 30 M), the remaining NS funnel curls and 
expands into a long, low-density tail while accretion con- 
tinues. During the tail phase, the accretion rate slows 
until « 99% of the NS rest mass is inside the BH (see 



Figs. 02 and O. 

In case B (q = 2), the BH is 33% smaller and X/"" 161 
reaches « 180° at t' ps 40M, at about the time when 
the densest part of the NS core falls into the BH. The re- 
maining material in the funnel begins to evolve into a tail 
around this time, when only w 55% of the NS rest mass 
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FIG. 2: Normalized violation of the ^-component of the mo- 
mentum constraint (top panel) and the Hamiltonian 
constraint \\TL\\ (bottom panel); see Eqs. (|40|l and (|41[) 



FIG. 4: Irreducible mass of the BH for sequence Al. Here M 
is the total ADM mass. 
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FIG. 3: Rest (baryon) mass fraction located inside the BH 
apparent horizon versus time for sequence A. Note that time 
is shifted {if = t — t-pc) so that the first contact occurs at 
t' = 0. The inset shows the rest mass fraction outside the BH 
horizon at late times. 



has fallen into the BH, and the accretion rate decreases 
considerably. Eventually the tail (Fig. [5]) is swallowed by 
the BH, but the accretion time scale is about twice as 
long as the q = 3 case (Fig. [5]). 

At t' ~ 30Af in case C (q = 1), a low-density re- 
gion of matter develops ahead of the higher-density fun- 
nel region as X/ nnel becomes larger than 180°. As the 
majority of the NS material passes through the high- 
density funnel, this low-density overshoot grows in size 
away from the horizon while wrapping quickly around 
the BH, and finally smashing into the high-density fun- 
nel region (Fig. [5]) before falling into the BH. Once the 
low-density overshoot has been accreted (t' w 125M), 
only about 15% of the NS rest mass remains outside the 
BH. The remaining NS matter is accreted at about the 
same rate as in the q — 3 case. 

We model our initial NS by a polytropic EOS P = 
kpq with r = 2. During the evolution, shocks develop 
and the matter heats up, resulting in an increase in the 
parameter K = P/(kpq) from its initial value of unity. 
Fig. [7] plots contours of K at different points in time for 
case C. Notice that K is clearly larger than unity in the 
low-density overshoot region (upper- left panel of Fig. [7]) , 
and where the overshoot region smashes into the higher- 
density funnel region (upper-right panel of Fig. [7J . We 
also observe shock heating as the final bit of NS matter 
is accreted (lower- left panel of Fig. [JJ. 

A variety of physical effects are not modeled in the 
post-shocked, semi-degenerate, nonequilibrium nuclear 
matter arising in these simulations. For example, trans- 
port due to photon and neutrino radiation is not mod- 
eled, so accurate measurements of temperature T are not 
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Case Al-hi 1 = 50. 8M Case B 1 = 59. 9M Case C 1 = 91. 9M 




X/M X/M X/M 



Case Al-hi t=132.9M Case B t=142.1M Case C t=173.9M 
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FIG. 5: Snapshots from runs Al-hi, B, and C, compared at first contact (upper panels) and at a moment in time At « 75M 
later (lower panels) . Contours are defined as in Fig. [1] 
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FIG. 6: Rest (baryon) mass located inside the BH apparent 
horizon versus t'/M for runs Al-hi, B, C, D. 



possible. Moreover, we are not employing a realistic hot 
nuclear EOS. However, we can very roughly estimate T 
from the specific energy density e. For a poly tropic equa- 
tion of state, the "cold" contribution e co id is 



Ccold = / icoldd(l//?o) = 



r - 1 



(55) 



where P C oid = ^Po ■ We now define the thermal contribu- 
tion to the specific energy density as eth = e — e coid and 
compute the thermal contribution according to 



1 P 



e cold 



J'-l 



r - 1 p a r - 1 



Po 



(K-l)e c 



>ld > 



(56) 



where we have used (|25[) to express e in terms of P and 
Po- 

To estimate T, we assume that we can model the tem- 
perature dependence of eth as 



eth 



3kT + _^aT i 
2m n po 



(57) 



(compare [92j), where m n is the mass of a nucleon, k 
is the Boltzmann constant and a is the radiation con- 
stant. The first term represents the approximate ther- 
mal energy of the nucleons, and the second term ac- 
counts for the thermal energy due to radiation. The 
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Case C t = O.OM Case C t = 91.9M 




-10-8 -6 -4 -2 2 4 6 8 10 -10-8 -6 -4 -2 2 4 6 8 10 



X / M X / M 

Case C t = 257.2M Case C t = 283.6M 




-10-8 -6 -4 -2 2 4 6 8 10 -10-8 -6 -4 -2 2 4 6 8 10 

X / M X / M 



FIG. 7: Snapshots of equatorial contours of the polytropic constant K at selected times for case C. Contours are spaced linearly, 
so that AK — 0.89 from K = 2 to K = 10. Values of K > 1 result from shock heating; K = 1 for adiabatic flow. 
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factor / reflects the number of species of ultrarelativis- 
tic particles that contribute to thermal radiation. When 
T <C 2m e /k ~ 10 10 K, where m e is the mass of elec- 
tron, thermal radiation is dominated by photons and 
/ = 1. When T ^> 2m e /fc, electrons and positrons be- 
come ultrarelativistic and also contribute to radiation, 
and / = 1 + 2 x (7/8) = 11/4. At sufficiently high tem- 
peratures and densities (T > 10 11 K, p > 10 12 g cm -3 ), 
neutrinos are generated copiously and become trapped, 
so, taking into account three flavors of neutrinos and anti- 
neutrinos, / = 11/4 + 6 x (7/8) = 8. For a heated region 
in the tidally disrupted NS, which possesses a typical 
density of 5 x 10 13 g cm -3 , and K ~ 7, we obtain a tem- 
perature of T ~ 2 x 10 11 K. Although these temperatures 
are sufficient to produce copious neutrino emission, our 
low disk mass would limit the overall neutrino energy to 
E v < 10 49 ergs (following the approximate scalings de- 
rived by |93| | from numerical models of BH disks). This 
limits the total 7-ray annihilation energy to < 10 48 ergs 
assuming 10% efficiency, which may not be sufficient to 
power a SGRB. We note, however, that SGRB produc- 
tion may not require a long-lived massive disk, since the 
actual emission mechanism remains poorly understood. 
Instead, we merely require that sufficient thermal energy 
be produced to power the burst itself. So until a more de- 
tailed model of SGRB generation from BHNS mergers is 
developed, any assessment of these simulations regarding 
SGRBs is tentative at best. 

In Fig. we plot snapshots of case D, the low- 
compaction, q — 2 mass ratio case. At the time of first 
contact, the accretion funnel is much narrower than in 
the high-compaction cases of Fig. [5l Further, unlike any 
of the high-compaction cases, ^ nnel surpasses 180° at 
t' w 165M - about 60M before the highest-density region 
of the NS has been accreted. After x ^™ ci = 180°, a case 
C-like overshoot develops (lower-left panel of Fig. [5]), but 
instead of smashing into the higher-density funnel and 
quickly falling in to the BH, the overshoot gently merges 
with the funnel to create a short-lived disk-like structure 
(lower-right panel of Fig. [5]) that later falls into the BH. 

The final disk masses we measure for each run are listed 
in Table IIII1 along with an estimate of the final (Kerr) 
BH spin. To calculate the latter quantity, we take the 
ratio of the polar to the equatorial circumference of the 
apparent horizon, C r = C p /C e and use Eq. (5.3) of [§4| 
to solve for the dimensionless spin a = a/M-u, where 

(M irr /a) ^2(1 



At- 



tn 

mass 



(1 - Vl - a 2 ) is the Kerr BH ADM 



C r = 



E - 



(58) 



Here E(x) is the complete elliptic integral of the second 
kind. In all cases, we find disk masses of less than 2.8% of 
the original NS mass. The final (Kerr) BH spin is roughly 
a/At = 0.5, 0.64, 0.8, and 0.5 for cases A, B, C, and D, 
respectively. Here we use our finding that At = 
to good approximation. The first of these agrees well 



with similar results of ST; this is the only case for which 
a meaningful comparison is possible, given the adopted 
mass ratios. 



C. Gravitational Wave emission 

The gravitational wavetrain from a compact binary 
system may be separated into three qualitatively different 
parts: the inspiral, merger, and ringdown. We describe 
each briefly before discussing our numerical results. 

During the inspiral phase, which takes up most of the 
binary's lifetime, GW emission circularizes the orbit and 
gradually reduces the binary separation. At the large 
binary separations during the inspiral stage, finite-size 
effects associated with the NS are unimportant, and post- 
Newtonian (PN) techniques are sufficient to describe the 
evolution. At present, the binary orbit dynamics is deter- 
mined to 3.5PN order (e.g. [95|) and the corresponding 
GW emission is computed to 2.5PN order [96|, [13] (but 
see also [Hj]). Even after finite-size corrections become 
relevant, quasi-equilibrium sequences allow for a deter- 
mination of the binding energy as a function of orbital 
frequency, from which the GW energy spectrum dE/df 
may be calculated, following the techniques described in 
[99| |. This method was used in [35| to determine the ap- 
proximate energy spectrum from the sequences we use as 
initial data in this work. 

Once the binary nears the ISCO, or the point where 
tidal disruption begins, the orbit decays rapidly and the 
GW emission changes character. In particular, devia- 
tions from point mass behavior typically result in a sharp 
decline in the energy spectrum. We note that these 
"break frequencies" marking the onset of instability sys- 
tematically occur at lower frequenci es fo r BHNS than 
for BHBH binaries (see discussion in 100]), especially in 
cases where tidal disruption occurs. This can be seen 
clearly in (36[, noting that the tidal disruption branches 
do not exist for BHBH systems. 

Finally, we expect a phase of quasinormal ringing of 
the BH, since it is distorted by the merger. This emission 
typically results in a higher frequency peak in the energy 
spectrum, with an amplitude determined by the total 
distortion induced on the BH by the merger. 

In Fig. [51 we plot the GW strains along the polar axis 
of the binary for case Al-hi, using both the ip^ formalism 
(solid curves; Eq. (|49j) ) and the Z-M formualtion (dashed 
curves; Eq. (|4"6")0 . In both cases the waveforms are ex- 
tracted on a sphere of physical radius r ex = 34. AM, and 
modes up to and including I = 4 are used in the cal- 
culation. We add suitable integration constants when 
computing the waveforms from both Z-M (odd-parity 
modes) and ipi formalisms to minimize offsets in the 
time-averaged h + and h x ■ We see that there is some dis- 
agreement at early times as spurious gravitational radia- 
tion present in the initial data propagate outward. Once 
this "junk" radiation has left the numerical grid, the two 
independent methods yield results that are in very good 



15 



12 



Case D 



>- 



4 


-4 E- 



12 



T 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



I I 



t = O.OM 



iii 



0.5C: 




I U U U LLL 



12 



Case D 



-4 



T 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



t=159.8M 



iii 



0.5C: 




I I 



U U U LL 



-12 

12 -8 -4 4 8 12 -12 -8 -4 4 8 12 

X / M X / M 




8 12 




8 12 



FIG. 8: Snapshots at selected times from case D, the low-compaction, q — 2 mass ratio case. The contours represent the 
density in the orbital plane, plotted logarithmically with four contours per decade, with greyscaling added for clarity. Arrows 
represent the velocity field in the orbital plane. The minimum contour value in each frame is ftpo(min)= 10~ 4 , or po(min) = 
1.6 x 1O 12 (1.4M /M o ) 2 g cm . The maximum initial NS density is npo = 0.058. We specify the black hole AH interior 
in each snapshot with a filled black circle. In cgs units, the total ADM mass for this case is M = 2 x 1O" 5 (M O /1.4M ) 
s= 6(M o /1.4M )km. 



agreement, even though they are calculated using differ- 
ent sets of metric components. 

During the inspiral phase, the GW frequency and am- 
plitude sweep upward until the point at which the NS 
begins to be disrupted by the BH, at t' = 0M. As accre- 
tion progresses , there is a gradual but steady downturn 



in the amplitude while the GW frequency continues to 
sweep upward. Finally, from t' = 50M onward, after the 
vast majority of the NS matter has been accreted, we see 
a very weak ringdown signal, at amplitudes significantly 
less than those seen in either BHBH ([13, 11 QjjlJ) or 
NSNS ([H) mergers. 



16 



TABLE III: Final results from each of our simulations. We list the fractional rest (baryon) mass outside the horizon / out = 
Mo(r > taw) /Mo and the dimensionless spin of the BH a/M [see Eq. ([58}] at the end of our simulation. Also shown are the 
radiated energy, angular momentum, and linear velocity "kick" resulting from GW emission, the former two normalized to the 
binary's initial total ADM mass and the latter in km/s. For the GW quantities, entries without (with) parentheses are derived 
from the Z-M (1/14) formalism. Note that in case D, the GW data are not accurate enough to obtain reliable estimate of A_Egw, 
AJgw and the kick velocity. 



Case 


font 


a/M AE G w/M 


AJ G w/M 2 


Kick velocity (km/s) 


Al-lo 


<l% 


w 0.52 0.60% (0.70%) 


5.2% (6.3%) 


21 (21) 


Al-med 


<l% 


ps 0.52 0.79% (0.75%) 


6.6% (6.2%) 


39 (20) 


Al-hi 


<l% 


ps 0.52 0.65% (0.77%) 


5.4% (6.4%) 


46 (19) 


Al-farbc 


<l% 


w 0.52 0.74% (0.74%) 


6.0% (6.1%) 


15 (16) 


A2 


<2% 


ps 0.52 0.72% (0.66%) 


7.3% (6.5%) 


49 (24) 


A3 


<2.8% 


ps 0.48 0.67% (0.86%) 


7.9% (5.1%) 


18 (30) 


B 


<1% 


ps 0.64 0.59% (0.53%) 


6.2% (6.2%) 


67 (48) 


C 


<1% 


ps 0.80 0.39% (0.32%) 


5.4% (4.7%) 


22 (25) 


D 


<1% 


ps 0.5 
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FIG. 9: Gravitational wave signal from case Al-hi, calculated 
using ip4 (solid line) and Z-M (dashed line). We show both 
polarizations as seen by an observer looking down the initial 
polar axis, h+ (top panel) and h x (bottom panel). The scale 
factor D is the distance to the binary. 



As discussed in Sec. lIV A} the numerical resolution has 
little effect on the mass accretion rate in sequence A. 
Correspondingly, we find that the tidal disruption sig- 
nature in sequence A waveforms is largely resolution- 
independent. For runs Al-hi, Al-med, and Al-low, we 
estimate the GW frequencies at t' — r ex = to be 
Mfl GW = 2vrAf/ G w = 0.188, 0.180 and 0.190, respec- 
tively. These frequencies are slightly higher than twice 
the orbital frequency value at tidal disruption found in 
[3^ . Mf2 or b ps 0.07, as is expected since first contact 



occurs slightly after the onset of tidal disruption. 

A similar pattern is observed in cases B and C (see 
Fig. [TO]). The ringdown amplitude grows relative to the 
overall signal strength as q is reduced from 3 (sequence A) 
to 1 (case C). However, in all cases the amplitude is sig- 
nificantly smaller than the comparable BHBH ringdown 
signal; we discuss this issue in more detail below. 

The low compaction NS in case D implies a larger NS 
radius, so a larger initial binary separation was required 
than for the other cases. Thus, case D required many 
more grid light-crossing times until merger and ringdown. 
As a result, late-time normalized Hamiltonian constraint 
violation \\H\\ increased to ps 8%, leading to an inac- 
curate late-time waveform. We therefore truncate the 
waveform after t' — r ex — 200M . Based on our analysis 
of sequence A and the known high computational cost 
of case D, accurate simulations of case D would require 
higher resolution and more distant outer boundaries than 
is practical, given our computational resources. Though 
we can still evolve the matter and trajectory of the BH re- 
liably, we find that since the waveforms are manifested as 
small perturbations on the background spacetime, they 
are greatly affected by constraint violations. 

We tabulate the GW energy loss Ai?cw and angular 
momentum loss AJqWj as well as the measured kick ve- 
locity imparted to the BH in the rightmost columns of 
Table Mil Quantities without parentheses are derived 
from the Z-M formalism waveforms, and those in paren- 
theses are derived from our -0 4 -based waveforms. We 
compute AEqw, AJqwi and the kick velocity at 5 radii 
in the range ps 31 — 37M for all cases except Al-farbc 
(where the radii span ps 47 — 76M). The values shown 
in Table IIIII are obtained by Richardson extrapolation of 
the data to r — > 00. In general, we find good agreement 
between the two GW measurement methods, especially 
for the case with more distant outer boundaries, which 
satisfied the constraints best at late times. Based on the 
variations of results with different resolutions (for cases 
Al), different GW extraction radii and in the two GW ex- 
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FIG. 10: Gravitational wave signal from cases B, C, and D, 
calculated using the Z-M formalism. We show both polariza- 
tions, h+ (solid) and h x (dashed) along the polar axis of the 
binary. The scale factor D is the distance to the binary. 



traction methods, we estimate that our tabulated A_Bqw 
and AJ is accurate to about 20%, whereas the error in 
kick velocity may be as much as 50%. In case D, the GW 
data are not accurate enough to provide reliable data for 
energy, angular momentum losses and kick velocity. 

Compared to previous simulat ions of merging BHBH 
systems with the same mass ratio jl02| | , we find that while 
the radiated energies, angular momenta and kick veloci- 
ties are significantly lower in our runs because tidal dis- 
ruption suppresses the GW signal, the final BH spins are 
comparable, within our uncertainties, to BHBH values. 

In a previous work [3c| . we showed that our -04 mea- 
surements converged to second-order with numerical res- 
olution. In Fig. 111! we perform a similar demonstra- 
tion, but with the Z-M formalism. In this figure, we 
plot the real component of the ^H en mode for cases 
Al-hi, Al-mcd, and Al-low, which differ only in the 
numerical grid spacing. The three waveforms are plot- 
ted in the top panel and show good agreement. No- 
tice that the waveform amplitudes are weakly depen- 
dent on the numerical resolution, but only at the level 
of a few percent for our higher resolution runs. In 
the bottom panel, we show differences between pairs 
of runs, rescaling the higher-resolution case by a nu- 
merical factor that assumes second-order convergence, 
(166~ 2 - 220~ 2 )/(220- 2 - 270~ 2 ) = 2.25, finding agree- 
ment. Although our spatial differencing scheme for the 
fields is fourth-order accurate, our HRSC scheme in un- 
shocked regions is only second-order accurate. This, 
along with the appearance of shocks (which are only 



FIG. 11: Numerical convergence of gravitational wave sig- 
nals for cases Al-hi, Al-med, and Al-lo. In the top panel, 
we show the I — m — 2 component of the even-parity mode 
^eSen f° r the three waveforms, noting that while they remain 
in phase with each other, we see overall amplitude differences 
on the order of several percent. In the bottom panel, pair- 
wise differences between the waveforms are plotted, with the 
higher-resolution pair rescaled to demonstrate second-order 
convergence. 

first-order convergent) limits the convergence order of our 
waveforms over time. As a result, while we can perform 
fourth-order time integrations with our code, we gener- 
ally prefer second-order time integration since it is faster 
and does not result in a significant loss of accuracy. 

Although the I = 2, m — 2 mode is the dominant 
component of the radiation, we measure all spin- weighted 
spherical harmonic components up to and including I = 
4. In Fig. [12l we show the mode decomposition of ip4 as 
a function of time for all non-negligible contributors. In 
the top panel of the figure, components satisfying / = m 
are plotted, including the dominant I = m = 2 mode. 
Notice that the modes satisfying I = m + 1 in the middle 
panel possess an amplitude that is at most 15% of the 
total strain at any given moment. For completeness, we 
plot some of the other significant modes in the bottom 
panel, noting that while they are present in the initial 
passage of "junk" radiation, they play little or no role at 
later times. 

To study the detectability and qualitative features 
of our computed gravitational wavetrains, we calculate 
the effective GW wave strain in frequency space using 
Eq. (|5~4|) . We find that if an FFT is performed on these 
GW signals without modification, the initial burst of spu- 
rious junk radiation contributes significantly to the sig- 
nal, and the finite initial amplitude introduces a strong 
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FIG. 12: Mode decompostition of the real component of ip4 
for case Al-hi. In the top panel, we show the dominant I = m 
modes. In the middle panel, we plot the modes that satisfy 
I = m + 1, and in the bottom panel we show several other 
modes. 



aliasing signal across the entire frequency domain of in- 
terest. To fix these problems, we perform "surgery" 
between our numerical GW signal, at a point in time 
where the initial junk has passed through the GW ex- 
traction surface, and a post-Newtonian signal with the 
same q. We generate the restricted P N wa veform fol- 
lowing the same techniques as in, e.g., 100]. Defining 
v = (ttM/gw) 1 ^ 3 , w e take as an initial condition the 
value of v computed from the phase evolution of our nu- 
merical waveform and evolve backwards in time the fol- 
lowing set of equations, 



dip 
~dt 
dv 
It 



2v 3 
~M 



F(v) 



M dE{v)/dv ' 



(59) 
(60) 



where the binding energy per unit mass E(v) and radia- 
tion flux F(v) are taken from the PN calculations of [95| . 
We modify the amplitude of the PN signal to minimize 
aliasing, but the relative correction is in all cases less 
than 2%. Mismatches in the amplitude, frequency, and 
frequency sweep rate appear as oscillations in the energy 
spectrum near the surgery frequency. As a result, any 
peaks and troughs appearing at comparable frequencies 
in the energy spectrum should be viewed with skepti- 
cism unless demonstrated to be robust with respect to 
the surgery procedure. At both the beginning and end 
of the combined waveform, we add exponential damping 
terms to reduce aliasing, but this operation does not add 



power at frequencies of interest. 

The result of this operation for run A3 is shown in 
Fig. 1131 where units are set by assuming a NS rest mass 
of 1.4M Q . We plot the effective strain computed from 
both our combined waveform and our numerical signal 
alone. Notice that the signal follows the point-mass 
power law behavior up to frequencies of approximately 
/gw = 600 — 800 Hz, at which point disruption of the 
NS and its subsequent accretion dominates the signal. 
At frequencies above /gw = 1 kHz, however, there is ex- 
tremely little power in the waveform since the ringdown 
signal is so weak. These results are consistent with what 
we expect from the quasi-equilibrium configurations of 
[3^ | . who find that /gw = ^orb/ 7 !" ~ 800 Hz for the con- 
figuration in question (A//gw — 0.022 in dimensionless 
units). They also agree roughly with results of [2l|, who 
find a similar pattern of steep decline above the tidal 
disruption value. 

Each BHNS merger spectrum is compared to a BBH 
merger spectrum, taken from Eqs. (4.12)-(4.19) of [lOOj ]. 
noting that /i e ff(/) oc fA(f) in their notation. The com- 
parison is performed using a binary with the same masses 
as our BHNS case. Both of these curves lie above the 
advanced LIGO sensitivity band /iligo(/) = \J fSh(f), 
which we have taken from [103J. This result assumes 
a distance to either source of D = 100 Mpc, the dis- 
tance required to reach one merger per year assuming an 
overall rate of 10 mergers per megayear per Milky Wav- 
equiv alent galaxy (and a density of these of 0.1 gal/Mpc ) 
jl04l | . This distance is roughly that of the Coma cluster, 
and approximately five times the distance to the Virgo 
cluster. The difference in wave signal between BHBH 
and BHNS mergers is present in the advanced LIGO fre- 
quency band, but only marginally. It is clear that for 
more significant measurements of the difference between 
BHNS and BHBH inspirals and mergers, it would be ad- 
vantageous to make use of narrow-band detection tech- 
niques with advanced detectors. 

Fig. [TJ] contains plots of the GW spectrum for cases 
B, C, and D. As the mass ratio and NS compaction is 
varied, we see the expected differences in the apparent 
"break frequency" marking tidal disruption. This fre- 
quency may be estimated by M /break ~ (M/dtid) 3 ^ 2 ~ 
C 3/2 (l + g) v /(l + q)/q/-K, where Eq. (fTJ) has been used for 
dtid- This formula is consistent with Eq. (25) of [36| . As 
the value of q is lowered from 3.0 to 1.0, the break fre- 
quency rises. In addition, when we lower the compaction, 
we see a large decrease in the break frequency. These re- 
sults, which agree well with the empirical scalings derived 
in [361 ] , lend credence to the idea that if the individual 
masses of the binary components can be derived from the 
inspiral waveforms, the GW break frequency should pro- 
vide a relatively sensitive measurement of the NS radius. 
When combined with observations of the lower-frequency 
spectrum, these inferences may provide additional con- 
straints on the NS structure, including limits placed on 
the tidal Love number &2 (following the techniques de- 
scribed in [l05| ] for NSNS mergers). 
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FIG. 13: Gravitational wave spectrum for the case A3 BHNS 
merger compared to a BHBH merger with the same masses. 
The solid curve shows the combined waveform found by at- 
taching the restricted PN waveform to our numerical signal, 
while the dotted curve shows the contribution from the latter 
only, demonstrating the expected aliasing behavior resulting 
from FFTs of discontinuous functions. The dashed curve is 
the analytic fit derived by [lOOfl from analysis of multi-orbit 
BHBH inspirals, which maintain significantly more power 
at higher frequencies. The heavy solid curve is the effec- 
tive strain of the advanced LIGO detector, defined such that 

hhiGo(f) = y fSh(f )- To set physical units, we assume a NS 
rest mass of M = 1.4M . 



V. DISCUSSION 

In this paper we present our first fully self-consistent, 
dynamical simulations of relativistic BHNS binaries. We 
use results and evolution techniques that we have pre- 
viously developed in preparation for these simulations, 
including the initial data of [3a|. the "filling" of the 
black holes in these initial data [38| and the treatment 
of relativistic hydrodynamics in the context of the mov- 
ing puncture method [371 ] . Here we focus on irrotational 
BHNS binaries with mass ratios between q = 1 and 3. 

For the cases studied here, we find no more than ~ 
3% of the original NS matter remaining outside the BH 
at the end of the simulations. Such small disk masses 
lend support to the semi-analytic arguments presented 
in 53] , which suggested that virtually the entire NS will 
be accreted promptly by the BH. 

The simulations of ST, on the other hand, suggest 
larger disk masses than ours. The reason for this discrep- 
ancy remains unclear, but we suspect it may be caused 
by different initial data - both choices of parameters 
and/or computational approach. We note that the mass 
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FIG. 14: Gravitational wave spectra for the cases B, C and D 
BHNS mergers compared to BHBH merger. Conventions are 
as in Fig. [T3] 



accretion process depends on the initial binary separa- 
tion (compare runs Al, A2 and A3 in Fig. [3]), suggesting 
that the disk mass depends rather sensitively on the de- 
tails of the initial data. However, we cannot rule out that 
this dependence is instead a numerical artifact, caused by 
numerical errors that are growing over time due to outer 
boundaries that are too close to the strong-field region 
and other effects. Since we model NSs as T = 2 poly- 
tropes, and the disk mass is likely to depend sensitively 
on the NS EOS, firm conclusions about BHNS mergers 
as SGRB progenitors remain uncertain. 

Our next series of BHNS simulations will involve spin- 
ning BHs. Most formation scenarios for BHNSs favor 
spinning BHs, especially systems in which the BH spin 
and orbital angular momenta axes are nearly aligned 
[l06j |. Since the ISCO for a prograde BH lies at a smaller 
radius than that of a non-spinning BH, we expect that 
the tidal disruption of the NS around these BHs occurs 
farther from the ISCO. Thus, spinning BHs would likely 
lead to a more massive disk, but the magnitude of the 
effect and the scaling with respect to spin will need to be 
determined via numerical calculations (as suggested in 
[Hlj . for the q ~ 10 cases). Such calculations will enable 
us to probe in depth which areas in phase space are likely 
to serve as progenitors for SGRBs. 

We find that the GW signal resulting from our BHNS 
coalescences is attenuated at frequencies roughly equal to 
double the orbital frequency at which tidal disruption be- 
gins, as one would expect, confirming the fits described 
in [3Q|. The deviation between BHNS and BHBH in- 
spiral is visible in the advanced LIGO band for systems 
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with mass ratios q = 3 out to distances > 100 Mpc, 
within which volume some population syn thesi s calcula- 
tions predict ~ 1 BHNS merger per year |l04l |. Should 
the chirp mass determination, combined with higher or- 
der PN waveform phase effects, allow for an independent 
determination of the component masses of the binary, 
observation of the BHNS merger break frequency should 
give a solid estimate of the NS radius. Such effects are 
independent of the discussion of disk formation, since the 
GW signal is strongly suppressed after the onset of tidal 
disruption. 

We have performed a series of calculations for the q = 3 
mass ratio case (sequence A), in which we vary only 
grid parameters to determine the numerical resolution 
requirements for these BHNS mergers. Gross features of 
the hydrodynamics, such as the accretion rate onto the 
BH, seem insensitive to the numerical resolution on the 
grid, at least for this case, where the NS accretes fairly 
promptly. The frequencies of the waveforms at critical 
moments are similarly insensitive to resolution. Wave- 
form amplitudes, on the other hand, vary by a larger 
amount with respect to resolution, and are seen to be 
accurate only at times when the constraint violations re- 
main small. At late times, when they are largest, con- 
straint violations are dominated by finite boundary ef- 
fects, which can be greatly reduced by enlarging the phys- 
ical extent of the grid. Given the computational resource 
requirements for these simulations, it may be extremely 
costly to calculate waveforms accurate to a few percent 
using fisheye grids or similar fixed-mesh refinement tech- 
niques. 



We expect that high accuracy calculations spanning 
~ 10 orbits, as are currently performed in BHBH merg- 
ers, will require us to use adaptive mesh refinement 
(AMR) techniques. Our current technique includes a sin- 
gle high resolution grid that encompasses both the BH, 
NS, and surrounding strong-field region. Outside of this 
region is a transition to a lower resolution grid domain 
that extends to the outer grid boundary. With AMR, 
we will be able to focus this high resolution entirely on 
the two regions immediately surrounding the BH and NS. 
With significantly fewer gridpoints in the strong-field re- 
gion of the grid, we will be able to place more gridpoints 
in the low-resolution, weak-field region, thus extending 
our outer boundaries. In many ways, a relativistic hy- 
drodynamics code with AMR will likely become a key 
tool for simulating BHNS spacetimes. 
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